Abstract

Aim

Plants that host root‐symbiotic nitrogen‐fixing bacteria have an important role in driving terrestrial ecosystem processes, but N‐fixing ability is unequally distributed among plant taxa and ecosystems. Here we explore the large‐scale distribution of N‐fixing plant species worldwide.

Location

Global.

Time period

Present.

Major taxa studied

Vascular plants.

Methods

We estimated root‐symbiotic N‐fixing plant species diversity (as Shannon entropy) and relative richness (log‐ratio of N‐fixing to non‐fixing plant species) for c. 7,800 km2 hexagonal grid cells using the NodDB and Global Biodiversity Information Facility (GBIF) databases. Additionally, we explored the distributions of plant species associated with rhizobia, actinobacteria or cyanobacteria (relative to other plant species), and the relative richness of N‐fixing trees (log‐ratio of N‐fixing to non‐fixing trees). We related N‐fixing plant species distribution to environmental (climate, soil) and biogeographical (biome, realm) variables using multiple linear regression.

Results

N‐fixing plant diversity and relative richness showed unimodal relationships with latitude. Diversity of N‐fixing plants was highest in warm and wet climates, but in dry biomes and in Australasia. The relative richness of N‐fixing plants was highest in warm and dry climates, in tropical and temperate grasslands and in Eurasia. Plants associated with cyanobacteria were more widely distributed near the equator, while those associated with rhizobia were more prevalent at the edge of the tropics, and those associated with actinobacteria at higher latitudes (especially in boreal forests). The relative richness of N‐fixing tree species was highest in cold and dry areas and in boreal forests, with contrasting peaks in the Northern and Southern Hemispheres.

Main conclusions

The distribution of N‐fixing plant species exhibits regional hotspots and coldspots related to both environmental conditions and biogeographical history. Global N‐fixing plant distributions are different for the key root‐symbiotic bacterial groups. Information about N‐fixing plant distribution can improve global models of ecosystem functions and contribute to understanding how plants respond to global change.

1 INTRODUCTION

Symbiotic nitrogen (N) fixation reflects the ability of some plant species to fix gaseous N2 from the air with the help of (a) bacteria collectively referred to as rhizobia (e.g. Rhizobiaceae, Burkholderiaceae), (b) actinobacteria in the genus Frankia, or (c) cyanobacteria in the family Nostocaceae (Kuyper & De Goede, 2013; Vitousek et al., 2013). Although other, non‐symbiotic N fixation processes are globally more widespread (Crews, 1999), symbiotic fixation generally accounts for the largest proportion of N fixation per unit area (Boring et al., 1988). Symbiotic N‐fixing plants (N‐fixers) are often major ecosystem drivers (Munzbergova & Ward, 2002; Vitousek et al., 1987) due to their ability to enrich soil with N and by shaping various biotic interactions above‐ and belowground. There are global estimates of potential annual symbiotic N fixation and its distribution across different biomes (Amundson et al., 2003; Cleveland et al., 1999), as well as a recent global map for the distribution of symbiotically N‐fixing trees across most forest ecosystems (Steidinger et al., 2019). However, in order to understand the factors influencing primary productivity, driving ecosystem functions or estimating global nutrient cycles in terrestrial habitats, it is essential to have an overview of the global distribution of all vascular plants capable of symbiotic N fixation.

The current empirical information about the large‐scale distribution of N‐fixing plant species and the factors influencing the distribution patterns is fragmentary, as most published data concern N‐fixing trees. Two recent studies (Menge et al., 2019; Steidinger et al., 2019) mapped the distribution of the relative cover of N‐fixing trees and found that N‐fixing trees dominate in the tropics. Climate and soil conditions have been the main explanatory factors for N‐fixing tree distribution and latitudinal gradient. Globally, N‐fixing trees dominate in warmer areas and alkaline soils (Steidinger et al., 2019). Menge et al. (2019) also found that N‐fixing trees in the Americas are more prevalent in areas with high mean annual temperature and precipitation. However, large‐scale studies of tropical forests and savannas in Africa and South America (Pellegrini et al., 2016) as well as of Neotropical forests (Gei et al., 2018) have shown an increase of the share of N‐fixing trees in vegetation with decreasing precipitation.

Besides climate, it has been suggested that the evolutionary history of symbiotic N fixation may have shaped the current distribution of N‐fixing plants (Crews, 1999), although recent evidence for this is limited (Menge, Batterman, Hedin, et al., 2017; Menge & Crews, 2016). There are several biogeographical differences in N‐fixing plant distribution that cannot be explained by climate conditions alone. Menge et al. (2019) found a latitudinal gradient in N‐fixing tree distribution in the Americas, but no such trend was present in Asia despite similar climatic gradients. Because symbiotic nitrogen fixation is an expensive process requiring both large carbon investment for construction of bacteria‐harbouring nodules and energy for nitrogen reduction and bacterial life, N‐fixing plants typically grow in open habitats (Vitousek & Field, 1999). Thus, they are characteristic ecosystem components in early successional or disturbed habitats (Levy‐Varon et al., 2019; Vitousek et al., 1987).

N‐fixing trees are more common in tropical forests, where N bioavailability is high, but scarce in temperate regions that are expected to be more N‐limited and where N‐fixing plants should have a competitive advantage (Hedin et al., 2009; Vitousek & Howarth, 1991; Vitousek et al., 2013). This paradox is often explained by differential N fixation strategies in tropical and temperate regions—rhizobial N‐fixing trees (which are typically facultative, regulating N fixation to meet nutritional demand) dominate equatorward of 35°N, whereas actinorhizal N‐fixing trees (typically obligate, maintaining fixation regardless of soil fertility) dominate to the north (Menge et al., 2014; Sheffer et al., 2015). However, the paradox could also be explained by other growth forms rather than trees being more important in fixing N in temperate regions (Cleveland et al., 1999; Sprent et al., 2017).

Most large‐scale studies have used data that allow one to calculate or estimate the relative basal area of N‐fixing trees, but comparable results have also been reported for the relative genus or species richness of N‐fixing trees (e.g. Menge et al., 2019), indicating that abundance and diversity patterns are likely correlated for trees. Studies from grassland also demonstrate that the species richness of N‐fixing plants is a reasonable surrogate for N‐fixer abundance (Jin et al., 2013; Latz et al., 2012).

Global change alters ecosystems and can have an impact on the N cycle either via increasing N input for plants, changing microbial activity, or by changing the distribution and richness of N‐fixing plants in all terrestrial ecosystems (Hungate et al., 2003; Yuan & Chen, 2015). Therefore, knowledge about the global distribution of all N‐fixing plants is essential. We aim to address the global distribution of all root‐symbiotic N‐fixing vascular plant species using the best available data sources—the Global Biodiversity Information Facility (GBIF) and a recently published database of plant genera with root‐symbiotic nitrogen fixation (NodDB, Tedersoo et al., 2018) to quantify and map N‐fixing plant species richness and diversity globally. Specifically, (a) we explore latitudinal gradients in N‐fixer distribution and determine the main biogeographical and environmental drivers shaping these patterns; (b) we compare whether N‐fixing plants associated with different bacterial groups (rhizobia, actinobacteria, cyanobacteria) show different trends; (c) we test if analyses only considering tree species are representative of all vascular plants.

2 METHODS

2.1 Species data

We used the GBIF for vascular species distribution data. GBIF collates information about species and their geographical coordinates based on individual sites. We acknowledge that GBIF data may have spatial and taxonomic biases, but trust that it is the best available solution to be used as a surrogate for large‐scale species distribution at a global scale. We obtained 156,200,298 georeferenced vascular plant occurrences for 289,295 species (GBIF.org, 10 November 2017, GBIF Occurrence Download https://doi.org/10.15468/dl.4nqoev). We used the ‘CoordinateCleaner’ package (Zizka et al., 2019) in the R environment (R Core Team, 2019) to remove occurrences with potential geographical errors (i.e. data entry errors, or coordinates recorded for capitals, country centroids and biodiversity institutions). The cleaned GBIF data were then matched to the latest GBIF backbone taxonomy (GBIF Secretariat, 2019) using the ‘taxize’ package (Chamberlain & Szöcs, 2013; Chamberlain et al., 2019) in R to resolve taxonomy issues. Finally, our dataset included 266,074 species with GBIF occurrence data.

We used a recently published global database of plants with root‐symbiotic N fixation ability (NodDB, Tedersoo et al., 2018) that compiles all available information about plant genera associated with N‐fixing root‐nodulating bacteria and supplements this with phylogenetic information. The resulting database describes plant genera that are associated with rhizobia, actinobacteria (Frankia) or cyanobacteria (root‐nodulating Nostocaceae), and also includes genera that are likely nodulated based on their phylogenetic information. In our study, we included all the aforementioned genera as N‐fixers, the rest were assigned as non‐fixers (as suggested by Tedersoo et al., 2018). Note that since we considered only root‐symbiotic N fixation ability, plant species from the genus Gunnera that host symbiotic cyanobacteria in the lower parts of leaf petioles (Franche et al., 2009) are considered non‐fixers in our dataset. The NodDB is collated at the plant genus level, but since nitrogen‐fixing ability is relatively constrained within genera (Sprent et al., 2017), we assigned every species in our GBIF dataset a nitrogen‐fixing ability based on its genus. Prior to assigning species with nitrogen‐fixing ability, the taxonomy of NodDB was also matched to the latest GBIF backbone taxonomy (GBIF Secretariat, 2019) using the ‘taxize’ package (Chamberlain & Szöcs, 2013; Chamberlain et al., 2019) in R to account for nomenclature differences between datasets. Finally, out of 266,074 species with GBIF occurrence data, 17,386 were considered N‐fixers according to NodDB.

The diversity and relative richness of N‐fixers was estimated for discrete global grids based on equal area polygons (Icosahedral Snyder Equal Area Aperture 3 Hexagonal Grid, Sahr et al., 2003). We used a grid that divides the globe into c. 65,000 hexagons, each c. 7,800 km2 and whose centres are separated by c. 100 km (loosely resembling the scale of 1 degree latitude × 1 degree longitude cells). We used the ‘dggridR’ package (Barnes, 2018) in R for creating equal‐area grid cells. For the analyses, we only considered terrestrial grid cells where the total number (including N‐fixers and non‐fixers) of plant species records exceeded 100. The final dataset covers 8,197 grid cells.

For each grid cell, we estimated N‐fixing plant species diversity as Shannon entropy with the number of records within a grid cell as a surrogate for abundance measure following Chao et al. (2013) using the ‘iNEXT’ R package (Chao et al., 2014; Hsieh et al., 2019). This estimation of Shannon entropy accounts for differences in sampling intensity between grid cells by extrapolating diversity to an asymptote and obtaining the expected diversity at full sample coverage. We also quantified the relative richness of N‐fixing plant species in each grid cell as log‐ratio: ln(no. of N‐fixing species/no. of non‐fixing species). Similarly, we also considered the relative richness of different N‐fixer groups (log‐ratio of plant species associated with either rhizobia, actinobacteria or cyanobacteria relative to all other species). For rhizobia, we included plant species where the N‐fixing estimation in NodDB was ‘Rhizobia’ or ‘likely_Rhizobia’. When the associated bacterial type was missing in NodDB, the species were not included. In total, 16,534 plant species were associated with rhizobia, 387 with actinobacteria and 299 with root‐nodulating cyanobacteria (these included only species belonging to the order Cycadales in our dataset). Finally, we also considered the relative richness of N‐fixing tree species (log‐ratio of N‐fixing trees to non‐fixing trees) to explore N fixation among trees. Tree data were obtained from the GlobalTreeSearch (BGCI, 2019), listing 60,011 tree species (Beech et al., 2017). After harmonizing the tree data to the latest GBIF backbone taxonomy (GBIF Secretariat, 2019) using the ‘taxize’ package (Chamberlain & Szöcs, 2013; Chamberlain et al., 2019) in R and matching them with GBIF occurrence data, 48,155 tree species were included in our dataset, of which 3,439 were N‐fixers. For additional comparisons, we also quantified the relative records count of all N‐fixers (log‐ratio of the number of records of N‐fixing to non‐fixing plant species, indicating the relative frequency of N‐fixers) and raw N‐fixer species richness in each grid cell without considering sampling intensity.

2.2 Environmental data

We included environmental (climate and soil) factors to explore the conditions where N‐fixers are more prevalent. For each grid cell we compiled 19 current bioclimatic variables as well as elevation data from BioClim (Fick & Hijmans, 2017). In addition, soil total nitrogen (N), carbon and nitrogen ratio (C : N) and soil pH were extracted from the global‐scale WISE30sec database (Batjes, 2016). The database includes data for a range of soil depths, but we included data from the topmost layer only (0–20 cm). See Supporting Information Table S1 for the complete list of included environmental variables.

We performed principal component analysis (PCA) including all the climate and soil variables. Precipitation‐related variables (except precipitation seasonality) were ln‐transformed and all variables were standardized (mean = 0, SD = 1) prior to the PCA. We included the first four principal components in the subsequent analysis. Each of these components contributed > 5% to the total variation and together described c. 83% of the total variation of environmental conditions. The first PCA component was positively correlated with annual mean temperature, as well as with the minimum temperature during the coldest month and quarter and was negatively correlated with temperature seasonality (|r| > .90); the second component was correlated with precipitation during the driest month and quarter (> .80); the third component was most related to precipitation seasonality (r > .60); and the fourth component was negatively related to maximum temperature of the warmest month and mean temperature of the wettest quarter (|r| > .45). All soil variables (N, C : N, pH) were best related to the second PCA axis (r values .42, .50, −.64, respectively). See Supporting Information Table S1 for the complete correlation table and the proportion of variance described by each PCA axis.

2.3 Biogeographical data

We included biogeographical factors (biogeographical realm and biome) to account for regional and global evolutionary history in N‐fixing plant species distribution. We used the current World Wildlife Fund (WWF) terrestrial ecoregions to delineate each grid cell into a contemporary biogeographical realm and biome based on their central point (Olson et al., 2001). We combined ecologically similar biomes, resulting in seven biome combinations (Supporting Information Table S2). The biomes Mangroves (code 14), Lake (98), and Rock and Ice (99) were excluded due to lack of data and dissimilarity from other included biomes. The six realms included were Palaearctic (PA), Nearctic (NA), Afrotropic (AT), Neotropic (NT), Australasia (AA) and Indo‐Malaya (IM). Oceania (islands of the Pacific Ocean) and Antarctic were excluded due to lack of data.

2.4 Statistical analyses

All analyses were performed in the R environment (R Core Team, 2019). For the latitudinal gradients, we compared linear and polynomial fits to explore if N‐fixing plant distribution follows the absolute latitude linearly or has a peak at certain latitudes. We report the fit of the model with a lower Akaike information criterion (AIC) value. To explore how environment and biogeography predict N‐fixer distribution, we used the information theoretical approach (Burnham & Anderson, 2002). Using the ‘lm’ function in R, we included average grid cell elevation (ln‐transformed) and four environmental measures from the PCA as continuous variables, as well as biome and realm as categorical variables in multiple linear regression. All continuous variables were standardized (mean = 0, SD = 1) for better comparisons of the slope estimates. Separate models were fitted for the N‐fixer diversity, relative richness, raw richness and relative records count, as well as for the relative richness of rhizobia‐, actinobacteria‐ or cyanobacteria‐associated N‐fixers, and the relative richness of N‐fixing tree species. We considered all possible models and chose the best model for each response variable based on the lowest AIC value (Burnham & Anderson, 2002). If two or more models were considered equal (ΔAIC < 2), we still chose the model with the lowest AIC for representation of the results as the most parsimonious. We used the ‘MuMIn’ library (Bartoń, 2019) for the model selection, the ‘ANOVA’ function in the ‘car’ library (Fox & Weisberg, 2011) for estimating the F‐values for each variable in the best models, and the ‘emmeans’ library (Lenth, 2019) to compare the individual biome and realm values within models (Tukey’s test). Finally, we used variation partitioning (‘varPart’ function in the ‘modEvA’ library; Barbosa et al., 2016) to distinguish between the relative effects of environment (climate and soil conditions), biome or realm.

We also tested for correlations or associations between the included variables (Supporting Information Figure S1): Pearson’s r for continuous variables, partial association eta for continuous‐categorical variable combinations (‘etasq’ function in the ‘heplots’ R library; Fox et al., 2018), and Goodman and Kruskal tau for the two categorical variables (‘GKtau’ function in the ‘GoodmanKruskal’ R library; Pearson, 2016). None of the continuous variables used in the models showed collinearity (|r| < .4; Dormann et al., 2013). We did not include latitude in the models since it was highly correlated with PC1 (r = −.9).

3 RESULTS

3.1 Global patterns of N‐fixer plant diversity and relative richness

We found that the diversity (estimated as Shannon entropy) hotspots of symbiotic N‐fixing plants were concentrated on lower absolute latitudes, but high diversity was also found in areas with Mediterranean climate in Europe, Australia and the California‐Mexico region (Figure 1a,b). The relationship between absolute latitude and relative N‐fixing plant richness (log‐ratio of N‐fixing to non‐fixing plant species richness) was less obvious, with a high proportion of N‐fixers found around the 20 degree latitude, corresponding to the arid areas (Figure 1c,d). Relative frequency of N‐fixing plants (log‐ratio of records count of N‐fixers to non‐fixers) was correlated with N‐fixer relative richness (= .84) and both measures exhibited similar patterns (Supporting Information Figure S2a,b), whereas the raw N‐fixing plant species richness patterns were similar to the extrapolated Shannon diversity patterns (correlation r = .61; Supporting Information Figure S2c,d). The relative richness of N‐fixing trees to other trees was highest in the arid regions of Mexico and Australia (Figure 1e,f), and this was reflected in different latitudinal patterns in the Northern (peak around 20°N) and Southern (peak around 30°S) Hemispheres.

image
Global patterns and latitudinal gradient of N‐fixing plant diversity (Shannon entropy; a, b), relative richness (log‐ratio) of N‐fixing plant to non‐fixing plant species (c, d) and relative richness of N‐fixing trees to non‐fixing trees (e, f). On the maps, the data include terrestrial equal area grid cells (hexagons, c. 7,800 km2) that had a total number of records >100 (8,197 grid cells). The colour scale for each map is standardized to use quantiles. In the latitudinal gradient figures, points (± SD) represent averages for 5‐degree latitudinal bands and the line (confidence interval in grey shading) represents an overall polynomial fit using raw data

N‐fixing plants associated with rhizobia were distributed globally within the sampled space, and the relative richness of N‐fixing plant species associated with rhizobia (log‐ratio of N‐fixers associated with rhizobia to all other species) followed a trend similar to the relative richness of all N‐fixers (Figure 2a,b). The largest proportion (96%) of all the N‐fixing plant species was associated with rhizobia. N‐fixing plant species associated with actinobacteria were mostly limited to higher latitudes with the highest relative richness (log‐ratio of N‐fixers associated with actinobacteria to all other species) in North America (Figure 2c,d). In contrast, only 649 grid cells, mostly concentrated on lower latitudes, had at least one plant species associated with cyanobacteria (Figure 2e,f).

image
Global patterns and latitudinal gradients of the relative richness (log‐ratio) of N‐fixing plant species associated with rhizobia (a, b), actinobacteria (c, d) or cyanobacteria (e, f) to other plant species. On the maps, the data include terrestrial equal area grid cells (hexagons, c. 7,800 km2) that had a total number of records >100 (8,197 grid cells). Dark grey cells represent areas where plant species with the specific N‐fixing association were not reported. The colour scale for each map is standardized to use quantiles. In the latitudinal gradient figures, points (± SD) represent averages for 5‐degree latitudinal bands and the line (confidence interval in grey shading) represents a polynomial fit using raw data. In (f), one cell with absolute latitude >35° was combined with the 35° latitudinal band in the Southern Hemisphere, and four cells with latitude >45° were combined with the 45° latitudinal band in the Northern Hemisphere

3.2 Relationships with environment, biome and realm

Environment, biome and realm explained 21% of the N‐fixer plant diversity and 40% of the relative richness of all N‐fixers (Table 1). For the different bacterial associations, the models best explained the relative richness of N‐fixing plants associated with actinobacteria (42%), followed by rhizobia (40%) and cyanobacteria (26%). When considering only tree species, the models explained 64% of the relative richness of N‐fixing trees.

TABLE 1. Model parameters, adjusted R2 values and the number of considered grid cells (n) of the best models according to Akaike information criterion (AIC) for total N‐fixer diversity (Shannon entropy), the relative richness (log‐ratio) of N‐fixing plant to non‐fixing plant species, the relative richness of N‐fixing trees to non‐fixing trees, as well as the relative richness of plant species associated with rhizobia, actinobacteria or cyanobacteria to all other plant species
All N‐fixers Trees Associated with rhizobia Associated with actinobacteria Associated with cyanobacteria
Diversity Rel. richness Rel. richness Rel. richness Rel. richness Rel. richness
F value p F value p F value p F value p F value p F value p
Elevation 166.48 *** X 27.50 *** X 4.64 * X
PC1 (temperature) 159.60 *** 222.96 *** 111.81 *** 282.47 *** 221.19 *** 63.00 ***
PC2 (precipitation, low soil pH) 41.32 *** 1,429.33 *** 1,012.24 *** 1,290.18 *** 14.29 *** 29.00 ***
PC3 (precipitation seasonality) 6.51 * X 6.18 * X 60.48 *** 14.83 ***
PC4 (low max. temperature) 30.30 *** 31.28 *** 5.55 * 28.26 *** 3.78 n.s. 12.04 ***
Biome 27.01 *** 37.85 *** 26.27 *** 45.61 *** 10.92 *** 8.13 ***
Realm 169.72 *** 53.57 *** 388.35 *** 70.99 *** 105.20 *** 21.60 ***
Adj. R2 .21 .40 .64 .40 .42 .26
n 8,006 8,005 4,355 7,926 3,594 649

Note

  • PC = principal component. The explanatory variables marked as X were not included in the best model.
  • *** p < .001;
  • ** p < .01;
  • * p < .05; n.s. = p ≥ .05.

The largest proportion of N‐fixer diversity was explained by the realm (Table 1, Figures 3 and 5). The diversity was highest in Australasia and lowest in Indo‐Malaya regions. N‐fixer diversity was highest in dry biomes and lowest in boreal forests. N‐fixer diversity increased with elevation, PC1 (high temperature) and PC2 (high precipitation), but decreased with both PC3 (high precipitation seasonality) and PC4 (low temperature during the warmest and wettest period). Similar results were also found for raw N‐fixer species richness without accounting for sampling intensity (Supporting Information Table S3, Figures S3 and S4).

image
The effects of elevation, climate and soil conditions [principal component analysis (PCA) axes], biome and realm for N‐fixing plant diversity (Shannon entropy), relative richness (log‐ratio) of N‐fixing plant to non‐fixing plant species, and relative richness of N‐fixing trees to non‐fixing trees using partial regression plots from the best model (Table 1). The box‐plot figures represent median and percentiles and letters show significant differences between groups according to Tukey’s test (p < .05). Biomes are simplified classifications of World Wildlife Fund (WWF) biomes: alp.g = alpine grasslands; bor.f = boreal forests; dry = dry biomes; tem.f = temperate forests; tem.g = temperate grasslands; tro.f = tropical forests; tro.g = tropical grasslands. Biogeographical realms are AA = Australasia; AT = Afrotropic; IM = Indo‐Malaya; NA = Nearctic; NT = Neotropic; PA = Palaearctic

The relative richness of total and rhizobia‐associated N‐fixer plants responded similarly to all factors and were best described by the combination of environmental conditions and biome (Table 1, Figures 3–5-3–5). Both response variables increased with PC1 (high temperature) but decreased with PC2 (high precipitation) and PC4 (low temperature during the warmest and wettest period). The log‐ratio of total N‐fixing plant species richness relative to non‐fixing plant species richness as well as the log‐ratio of rhizobia‐associated plant species richness relative to all other plant species richness were highest in temperate and tropical grasslands and lowest in boreal forests. Both the relative richness of total and rhizobia‐associated N‐fixing plant species were also highest in the Palaearctic and Indo‐Malaya realms. Again, the relative records count of N‐fixers followed similar trends to the relative richness of N‐fixers (Supporting Information Table S3, Figures S3 and S4), but was also positively affected by PC3 (precipitation seasonality). The relative records count of N‐fixers was also highest in the Palaearctic realm (as the relative richness), but lowest in the Indo‐Malaya realm.

Relative richness of N‐fixers associated with actinobacteria had the largest proportion of variation explained by environment, biome and realm (Table 1, Figures 4 and 5). The relative richness of actinobacteria‐associated N‐fixers decreased with PC1 (high temperature) as well as with PC2 (high precipitation), but increased with both elevation and PC3 (high precipitation seasonality). The relative richness of N‐fixers associated with actinobacteria was quite uniformly distributed between biomes (highest value in temperate forests) and it was equally high in the Australasia, Indo‐Malaya and Nearctic realms.

image
The effects of elevation, climate and soil conditions [principal component analysis (PCA) axes], biome and realm for the relative richness (log‐ratio) of plant species associated with rhizobia, actinobacteria or cyanobacteria to other plant species using partial regression plots from the best model (Table 1). The box‐plot figures represent median and percentiles and letters show significant differences between groups according to Tukey’s test (p < .05). Biomes are simplified classifications of World Wildlife Fund (WWF) biomes: alp.g = alpine grasslands; bor.f = boreal forests; dry = dry biomes; tem.f = temperate forests; tem.g = temperate grasslands; tro.f = tropical forests; tro.g = tropical grasslands. Biogeographical realms are AA = Australasia; AT = Afrotropic; IM = Indo‐Malaya; NA = Nearctic; NT = Neotropic; PA = Palaearctic

Relative richness of N‐fixer plants associated with cyanobacteria was best explained by the realm (Table 1, Figures 4 and 5). It increased with PC1 (high temperature) and decreased with PC2 (high precipitation), PC3 (high precipitation seasonality) and PC4 (low temperature during the warmest and wettest period). The relative richness of N‐fixers associated with cyanobacteria was quite uniformly distributed between biomes and realms (highest value in temperate grasslands and Indo‐Malaya realm), but entirely absent from the boreal forest biome.

image
Proportion of variation explained by climate + soil conditions, biome and realm according to the best model (Table 1) for N‐fixing plant diversity (Shannon entropy; a), the relative richness (log‐ratio) of N‐fixing plant to non‐fixing plant species (b) and the relative richness of N‐fixing trees to non‐fixing trees (c), as well as the relative richness of N‐fixing plants associated with rhizobia (d), actinobacteria (e) or cyanobacteria (f) to other plants

The ratio of N‐fixing tree species to non‐fixing trees was best explained by the realm (Table 1, Figures 3 and 5). It was negatively related to elevation, PC1 (high temperature), PC2 (high precipitation), PC3 (high precipitation seasonality) and PC4 (low temperature during the warmest and wettest period). The relative richness of N‐fixing trees was highest in boreal forests and in the Australasia and Neotropic realms.

4 DISCUSSION

Nitrogen addition experiments show that N limitation is globally widespread (LeBauer & Treseder, 2008). In order to better understand the factors driving ecosystem functions and nutrient cycles or limiting primary productivity, and to predict global change scenarios, it is important to estimate the global distribution of all N‐fixing plants. We used an extensive global dataset (> 150 million georeferenced records) to explore the global diversity patterns of vascular plant species capable of N fixation via root‐nodulating bacteria. Climate and soil data accompanied by simple biogeographical variables such as biome and biogeographical realm explained up to 64% of the variation in N‐fixing plant global distribution. While both N‐fixer diversity (estimated as Shannon entropy) and relative richness (log‐ratio of N‐fixing to non‐fixing plant species richness) showed a unimodal relationship with absolute latitude, the global distribution of these variables varied owing to different responses to environmental (climate and soil) and biogeographical factors (biome and realm). Remarkably, the relative richness of N‐fixers associated with different root‐symbiotic bacterial groups showed contrasting global distribution patterns.

The global distribution of the relative richness of N‐fixing trees (log‐ratio of N‐fixing to non‐fixing tree species) is in agreement with the recent estimates about most forest ecosystems (Steidinger et al., 2019)—N‐fixers have the highest proportion around the 20 degree latitude band in the Northern Hemisphere and around the 30 degree band in the Southern Hemisphere, corresponding to the drier biomes. Both datasets were similar in terms of the proportion of N‐fixing tree species (7%). Other large‐scale studies have also shown that N‐fixing trees are more abundant in more arid habitats (Gei et al., 2018; Liao et al., 2017; Pellegrini et al., 2016, but see Menge et al., 2019). However, we found that the relative richness of N‐fixing trees was highest in colder regions, boreal forests and the Australasia and Neotropic realms. These results differ from the high relative abundance of N‐fixing trees at higher temperatures, in North America and Africa and the low dominance of N‐fixing in the boreal region as reported in Steidinger et al. (2019). Some of the differences may be due to the variation in the global patterns of tree abundance and tree richness. Since trees account for a larger proportion of global vegetation biomass (Erb et al., 2018), the results for relative richness of N‐fixing trees might contribute more to the global estimates of nutrient cycles, or be more relevant for forest ecosystems. However, trees constitute a minority of the total plant species richness in the world (Beech et al., 2017) and the ecological patterns evident for trees might not apply globally to all plant species.

Total N‐fixing plant diversity and relative richness also showed a unimodal latitudinal gradient and generally higher values in tropical regions, and there were no differences between the hemispheres. The relative richness of N‐fixers associated with rhizobia (log‐ratio of N‐fixing plant species associated with rhizobia to all other plant species) followed the same pattern, supporting previous findings of the high legume diversity in the tropics (Cleveland et al., 1999; Oliveira‐Filho et al., 2013). On the other hand, N‐fixers associated with actinobacteria were distributed mainly at higher absolute latitudes. This supports the idea that the obligatory association with Frankia may be a way to cope with N‐limited soils in temperate regions (Lu & Hedin, 2019; Menge, Batterman, Liao, et al., 2017; Menge et al., 2014). Cleveland et al. (1999) predicted that in boreal forests, much of the N fixation occurs via free‐living bacteria or bacteria associated with lichens or bryophytes. However, trees can also fix N in boreal forests via association with actinobacteria. Similarly to the turnover in the major symbiotic guilds (high abundance of arbuscular mycorrhizal plants near the equator, followed by high abundance of N‐fixers at mid‐latitudes and ectomycorrhizal plants at higher latitudes; Soudzilovskaia et al., 2019; Steidinger et al., 2019; Toussaint et al., 2020), there also seems to be a turnover in the N‐fixer types, with cyanobacterial associations dominating in the tropics, followed by rhizobial associations peaking at the 20–30 degree latitudinal band and actinorhizal associations dominating at higher latitudes.

The N‐fixing plant species diversity was high in areas with higher annual temperature and precipitation. The optimum temperature for biological nitrogen fixing is c. 25 °C (Houlton et al., 2008), with nitrogenase activity declining strongly below 22 °C (Ceuterick et al., 1978). Overall this pattern corresponds to the global pattern of vascular plant diversity—high annual energy input with constant water supply is associated with high plant diversity (Kreft & Jetz, 2007). High N‐fixer diversity may also help to explain the global vascular plant diversity patterns, because high N‐fixing activity promotes diversity of other species in surrounding areas indirectly through soil fertilization (Levy‐Varon et al., 2019, Vitousek & Walker, 1989).

Relative richness of all N‐fixers and of those associated with rhizobia increased with temperature and decreased with increasing precipitation. Pellegrini et al. (2016) showed that aridity favours richness and abundance of N‐fixing trees in tropical savannas and tropical forests. N‐fixers may use water more efficiently, which would give them a competitive advantage in arid areas and help to explain their distribution (Pellegrini et al., 2016). It has been suggested that N‐fixing legumes are more prevalent in dry rather than wet conditions, because N fixation and reduced leaflet size enhance the legumes’ water‐use efficiency and drought tolerance (Adams et al., 2016; Gei et al., 2018). Arid conditions are also distinguished by the low availability of total N in the soil (Delgado‐Baquerizo et al., 2013) making biological N fixation in those habitats physiologically more vital.

Tropical areas have not only large soil N storage (Amundson et al., 2003; Hedin et al., 2009), but also a high relative richness of N‐fixing plants. We were able to relate N‐fixing plant distribution to soil pH, and to some extent to N content and C : N ratio. While all the relative richness measures increased with alkaline soils (as in Steidinger et al., 2019), species diversity of N‐fixing plants was higher in more acidic soils. Biological N fixation itself can also lead to soil acidification, especially in warm and wet regions (Lu et al., 2014).

N‐fixer diversity was highest in the dry biome and tropical forests and grasslands. Relative richness of total N‐fixers as well as of those associated with rhizobia peaked in both tropical and temperate grasslands. Tropical grasslands exhibit the highest rates of biological N fixation per unit area (Cleveland et al., 1999) and this is in agreement with the observed high proportion of N‐fixing species (especially those associated with rhizobia) in tropical grassland habitats. In addition, the evidence of greater N‐fixer richness in dry habitats and higher elevations, which tend to be relatively open habitats, indicates that light availability is important to overcome the high energy cost of nitrogen fixation (Vitousek & Field, 1999).

Besides environmental constraints and biome‐related differences in N‐fixer distribution, differences were observed among biogeographical realms, pointing to historical dispersal barriers. In particular, Australasia had the highest diversity of N‐fixers, possibly owing to the high endemic diversity of root‐nodulating legumes (Sprent et al., 2017). While the realm explained a low amount of variation in the relative richness of total and rhizobia‐associated N‐fixers, the relative richness values were highest in the Palaearctic and Indo‐Malaya realms. It has been suggested that some evolutionary constraints on the N‐fixing plant species distribution and diversity may exist (Crews, 1999; Menge, Batterman, Liao, et al., 2017; Menge & Crews, 2016), which can explain differences between the relative richness of N‐fixers in realms or continents (Pellegrini et al., 2016). It is also possible that N‐fixing ability has been lost in species that during evolutionary history encountered conditions where the costs of N fixation outweighed its benefits (Griesmann et al., 2018).

So far, very little information exists about the global distribution of N‐fixing vascular plants associated with cyanobacteria. In our dataset as well, the data on cyanobacteria‐associated N‐fixers (only Cycadales) were extremely limited and their distribution was mainly concentrated in the tropical region. Nostocaceae bacteria are widely distributed in terrestrial environments, but can also grow symbiotically in plant tissues other than roots (e.g. soil‐buried lower petiole parts in Gunnera) and associate with other non‐vascular taxa (algae, lichen‐forming fungi, bryophytes; Franche et al., 2009), which were not represented in our dataset.

Our results show global trends in N‐fixing vascular plant species distribution, but more research is needed in data‐poor regions and to adjust the global results to specific ecosystems or regions. We accept that the current datasets have geographical and sampling biases, which we have tried to alleviate with data cleaning (Zizka et al., 2019) and by focusing on the relative richness or extrapolated diversity values. Relative values are more resilient to sampling biases and the estimation of Shannon entropy following Chao et al. (2013) should diminish uncertainties in heterogeneous datasets. We trust that our data are robust for highlighting global patterns in N‐fixers, but as more data of species distributions and densities become available, regional and local patterns could be unravelled to understand the relative role of N‐fixing plants in communities.

Importantly, our results represent the potential for N fixation rather than the actual rate. There are many mechanisms that can inhibit N‐fixing plants from overcoming nitrogen limitation at local scales, including energetic constraints (Dovrat & Sheffer, 2019; Sheffer et al., 2015), co‐limiting nutrients (Yuan & Chen, 2015), ecological constraints (e.g. fire, preferential grazing; Cech et al., 2010), competing and pathogenic microorganisms (van der Heijden et al., 2016), and human impact (Vitousek et al., 2013). In many ecosystems, it is also important to consider N fixation by free‐living microbes, bacteria inhabiting other vascular plant organs beside roots (e.g. leaf petioles in Gunnera), or bacteria associated with bryophytes and lichens (Cleveland et al., 1999; Franche et al., 2009).

The distribution of N‐fixing plant species is not globally uniform, but the share of N‐fixers exhibits regional hotspots and coldspots. Both N‐fixer diversity and relative richness are highest at around 20 degrees latitude. Biome‐wise, the greatest role of total and rhizobia‐associated N‐fixers is in temperate and tropical grasslands rather than forests. These ecosystems might be the most N‐limited, although the relationship between the presence of N‐fixers and N deficit is not straightforward (Vitousek et al., 2010). Moreover, the importance of the different types of N‐fixers changes geographically—cyanobacterial associations being represented in the tropics, actinorhizal ones in the higher latitudes, and the relative richness of rhizobia‐associated N‐fixers being more prevalent at the edge of the tropics and decreasing towards the temperate zones. Although the current work gives only the first insight into the distribution of vascular plant N‐fixers and the underlying factors, it is obvious that the heterogeneous distribution of N‐fixers has to be accounted for when estimating terrestrial productivity, ecosystem functions and nutrient cycles, as well as in global change models.

ACKNOWLEDGMENTS

We thank Brody Sandel and Lars Hedin for their valuable comments. This study was funded by grants from the Estonian Research Council (IUT 20‐28, IUT 20‐29) and by the European Regional Development Fund (Centre of Excellence EcolChange).

    BIOSKETCH

    This research was done within the Centre of Excellence EcolChange. EcolChange (http://ecolchange.emu.ee/en/) is a collaborative project of two universities, the Estonian University of Life Sciences (one top research team) and University of Tartu (four top research teams and an individual researcher). The research of EcolChange covers the field of global change ecology from molecular to biome‐level responses.